Unveiling new perspective of phylogeography, genetic diversity, and population dynamics of Southeast Asian and Pacific chickens

The complex geographic and temporal origins of chicken domestication have attracted wide interest in molecular phylogeny and phylogeographic studies as they continue to be debated up to this day. In particular, the population dynamics and lineage-specific divergence time estimates of chickens in Southeast Asia (SEA) and the Pacific region are not well studied. Here, we analyzed 519 complete mitochondrial DNA control region sequences and identified 133 haplotypes with 70 variable sites. We documented 82.7% geographically unique haplotypes distributed across major haplogroups except for haplogroup C, suggesting high polymorphism among studied individuals. Mainland SEA (MSEA) chickens have higher overall genetic diversity than island SEA (ISEA) chickens. Phylogenetic trees and median-joining network revealed evidence of a new divergent matrilineage (i.e., haplogroup V) as a sister-clade of haplogroup C. The maximum clade credibility tree estimated the earlier coalescence age of ancestral D-lineage (i.e., sub-haplogroup D2) of continental chickens (3.7 kya; 95% HPD 1985–4835 years) while island populations diverged later at 2.1 kya (95% HPD 1467–2815 years). This evidence of earlier coalescence age of haplogroup D ancestral matriline exemplified dispersal patterns to the ISEA, and thereafter the island clade diversified as a distinct group.


Results
Haplotype variation and genetic diversity. We analyzed complete mtDNA control region sequences of chickens from Cambodia (n = 173), Laos (n = 63), Thailand (n = 25), Myanmar (n = 78), the Philippines (n = 6), and Fiji (n = 24) generated in this study and including previously published sequences from the Philippines (n = 129), Indonesia (n = 10), and Pacific (n = 11). A total of 133 haplotypes were identified, with 70 variable sites consisting of 10 singletons and 60 parsimony informative sites. We documented 82.7% geographically unique haplotypes, while 17.3% of haplotypes were shared transregionally across SEA, suggesting high polymorphism among the studied individuals. Island populations (i.e., Philippine and Pacific chickens) accounted for 28% of all unique haplotypes identified, while 72% were unique to continental populations. Summary of observed polymorphic sites and haplotype variations are presented in Supplementary Tables S1-S2.
The indices of genetic diversity for each geographic population are shown in Table 1. Undoubtedly, haplotypic and nucleotide diversity was very high in SEA chicken populations. The mainland Southeast Asia (MSEA) chickens had higher total haplotypic diversity (Hd = 0.963 ± 0.005) and nucleotide diversity (π = 0.00782 ± 0.00398) than the Island Southeast Asia (ISEA) chickens (Hd = 0.942 ± 0.009; π = 0.00466 ± 0.00249), although no major Table 1. Genetic diversity indices of Southeast Asian and Pacific chicken populations estimated using complete mtDNA D-loop sequences. a Included published sequences of Pacific chickens (n = 11). b Combined sequences representing MSEA (Cambodia, Laos, Thailand, and Myanmar chickens). c Combined sequences representing ISEA (Philippine chickens and database sequences of Indonesian chickens; n = 10). N number of sequences, Ht number of haplotypes, Hd haplotype (gene) diversity, π nucleotide diversity. * p < 0.05; ** p < 0.01. www.nature.com/scientificreports/ differences were observed. The highest value of Hd and π was found in Thai chickens (with 72% RJF population in our dataset), whereas the least was observed in Pacific chickens. These results should be taken with caution, given the Pacific chickens have a relatively small sample size, and the Thai chicken dataset was predominated with RJFs (18/25) (Supplementary Tables S1-S2). Thus, the genetic diversity is usually higher than that of the domestic chicken populations. Remarkably, the Thai chickens had a high number of haplotypes (Ht = 19) in 25 individuals examined, suggesting a diverse population in the region. Similarly, intraclade diversity indices indicated high haplotype and nucleotide diversity of haplogroup D than all other major haplogroups classified in SEA and the Pacific chickens (Supplementary Table S3).
Phylogeography and genetic affinities of continental and island SEA chickens. The sequences generated in the present study and the reference sequences that represent chicken mtDNA control region-based haplogroup nomenclatures were used to reconstruct the matrilineal phylogeny (Supplementary Tables S1, S4). Pioneering molecular phylogenetic studies based on mtDNA control region and mitogenomes revealed fourteen haplogroups (A-I and V-Z) of chicken worldwide 7,8,19,31 . Divergent haplogroups D and V showed enigmatic phylogeny resolution and previously claimed to have been distributed in ISEA and Thailand, respectively 7,31,36,37 .
In this study, model-based maximum likelihood and Bayesian phylogenetic analyses produced concordant topologies and comparable branch lengths of the tree ( Fig. 1a; Supplementary Figs. S2-S3). Major clades have strong SH-aLRT and UFBoot supports for the ML tree and significant posterior probability support for the Bayesian tree. Minor differences involved only some rearrangements of terminals for haplotypes: Hap_60, Hap_61, Hap_62, Hap_66, and Hap_122, as they clustered with haplogroup D1 (i.e., sub-haplogroup D1a) in the ML tree ( Supplementary Fig. S2), while grouped with haplogroup D2 in the Bayesian tree ( Supplementary Fig. S3).
Phylogenetic analyses grouped the MSEA chickens into major haplogroups A, B, D, E, and F, with the evidence of newfound haplogroup V as a sister-clade of haplogroup C ( Fig. 1a; Supplementary Fig. S2). Haplogroup V, classified by ancestral mutation motifs A281G, T355C, and C363T, was further subdivided into two sub-haplogroups (Fig. 1b). Here, we documented evidence of sub-haplogroup V2 (classified by unique mutational motifs: C228T, A237G, C391T) only identified in Cambodian and Laotian domestic chickens, and one haplotype of Thai RJFs (G.  Table S1). Evidence of sub-haplogroup V1 observed predominantly in Thai RJFs (Hap_107, Hap_108, and Hap_110) shared commonality to the reclassified haplogroup V of RJFs in Thailand and Cambodia 31 ( Supplementary Fig. S3). Interestingly, both model-based phylogenetic trees revealed ancestral lineage of haplogroup D2 from MSEA chickens, mostly observed in Cambodian chickens (38.7%) and some low frequency of Laotian (7.9%) and Thai chickens (8.0%), while remaining undetected in Myanmar chickens ( Fig. 1a; Supplementary Fig. S2; Supplementary Table S1). Haplogroups A and B have wide geographical distribution all over SEA, while haplogroup F was prevalent among Myanmar chickens (34.6%), with some low frequency detected in Thai chickens (12.0%).
Likewise, ISEA chickens (i.e., the Philippines and Pacific) have a shared genetic affinity for predominant haplogroup D1. Godinez et al. 37 previously characterized this island sub-group as the "Philippine-Pacific subclade. " This sub-clade is defined by five unique mutational motifs, A281G, C296T, T306C, A342G, and G686A, and includes diagnostic motifs from the downstream region of the complete mtDNA control region sequence ( Fig. 1b; Supplementary Table S1). These findings also correspond to the diagnostic motifs (SNPs: A281G, C296T, T306C, A342G) of Polynesian chicken ancient DNA 36 .
Consistent classification of the major mitochondrial lineages of SEA chickens was also depicted in the medianjoining (MJ) network analysis (Fig. 2). Notably, haplogroup V lineage was separated from haplogroup D and haplogroup C with nine and seven mutational sites, respectively. Within the haplogroup V lineage, newly classified sub-haplogroup V2 was separated from sub-haplogroup V1 with four mutational signatures. The geographicalspecific MJ network analyses exhibited a close transregional evolutionary relationship of MSEA chickens in major haplogroups except for haplogroup F, which was predominated in Myanmar chickens (Supplementary Fig. S4a-d). Similarly, the Philippine and Pacific chickens also shared closely related haplotypes classified under sub-haplogroup D1b ( Supplementary Fig. S4e-f).
Population structure and genetic differentiation. We carried out a multivariate approach to complement the phylogenetic analyses to assess further the relationships among and between geographical populations, including database sequences of East Asia, South Asia, and Middle East chickens (Supplementary Table S5). The result of the PCoA distinguished population substructure between mainland and island SEA chickens along the first two axes, which accounted for 52.09 or 52.53% variation ( Fig. 3a-b). A homogenous subgroup was observed within island populations, particularly among the Philippine and Pacific chickens (F ST = 0.06936), while MSEA populations showed more a diverse assemblage, consistent with the phylogenetic analyses and haplogroup variations. In addition, we documented closer relationships between Myanmar chickens and Yunnan chickens than any other Chinese chicken population. The pairwise F ST value confirmed that Myanmar and Yunnan chickens were not differentiated from each other (F ST = 0.00816; p value < 0.01). Meanwhile, within MSEA chickens, transregional population substructures were observed ranging from 0.06895 between Laos and Thailand to 0.19202 between Cambodia and Myanmar (Supplementary Table S6). Interestingly, Cambodian chickens were situated halfway between other continental populations and ISEA chickens, supporting the basal affiliations of identified ancestral matriline (i.e., sub-haplogroup D2) depicted in both ML and BI phylogenetic trees. The PCoA plot also indicates a significant genetic differentiation and substructure between East Asian chickens and South Asian-Middle Eastern chickens, ranging from 0.14938 to 0.77115 ( Fig. 3a; Supplementary Table S6). Similarly, we observed a close genetic affinity of Japanese and Korean chickens to the Chinese and MSEA chicken populations after removing South Asian and Middle Eastern chickens from the dataset (Fig. 3b). www.nature.com/scientificreports/ Hierarchical AMOVA revealed that the majority of the variations (i.e., 79.21% between ISEA and MSEA chickens and 79.74% between MSEA and EA chickens) could be attributed to within-population differentiation, specifically chickens distributed across Southeast and East Asia ( Table 2). Higher within-population variation was also observed within ISEA and Cambodian chickens. Likewise, no significant population genetic differentiation was found among groups of the island and mainland SEA chickens and groups of MSEA and East Asian chickens. These observed patterns of genetic differentiation from the partitioned variances among hierarchical groups reflect consistency established in the previous phylogenetic and PCoA analyses.
Demographic history and divergence time estimate. The simulations for neutrality tests indicated both MSEA and ISEA chickens deviated from neutrality (Table 1)  The maximum clade credibility (MCC) tree estimating the divergence time using a calibration method under an uncorrelated lognormal relaxed clock model revealed age estimates for biogeographically important nodes of haplogroups D and V in our dataset (Fig. 4c). The node age of macrohaplogroup CDV was estimated to be 6.67 kya with credibility intervals of 4235-7996 years (95% HPD). The coalescence age of sub-haplogroup D1b (PP = 1) was dated back to 2.1 kya (95% HPD 1467-2815 years) while diverging from the ancestral D-lineage approximately 3.7 kya (95% HPD 1985-4835 years). Haplogroup CV (PP = 0.80) diverged much earlier from macrohaplogroup CDV and coalesced around 5.5 kya (95% HPD 3116-7275 years) while succeeding divergence of haplogroup V (PP = 0.96) occurred around 3.9 kya (95% HPD 2125-5880 years). Newfound evidence of

Discussion
The timing and location of chicken domestication have been the subject of protracted debate worldwide and have stimulated several molecular studies using modern biological and zooarchaeological data 8,9,[11][12][13]20,31 . The consensus among researchers and several molecular studies confirmed that domestic chickens evolved from red junglefowl somewhere in South and Southeast Asia 7,8,11,28,31 , but identifying their exact geographic center of origin has been challenging 9,12,13 . Here, we present a comprehensive resolution of mitochondrial lineage diversity and phylogenetic analyses, population differentiation, demographic inference, and divergence time estimates of chickens in Southeast Asia and the Pacific region. Patterns of sequence variation indicated that chickens in the MSEA region have higher intrapopulation genetic diversity than island populations. The average genetic diversity   43 , and East African chickens (Hd = 0.638 ± 0.024; π = 0.00745 ± 0.00042) 32 . The substantial diversity of SEA chickens reflects the high matrilineal genetic variation documented in the major haplogroups, particularly haplogroup D with a large number of divergent haplotypes and haplogroup V, which has been detected only in Thailand, Cambodia, and Laos (Supplementary Tables S2, S3). However, we cannot invalidate the influence of RJFs samples on the overall genetic diversity as they reflect ancestral genetic variations. Divergent sub-haplogroups that retained ancestral variations were also observed in these lineages, likely due to geographic proximity to the center of domestication. These defined indices of biodiversity offer great opportunities for developing genetic improvement strategies, trait selection, effective management of genetic resources, and future conservation efforts 1,4,44,45 .
Pioneering molecular studies and DNA sources based on the hypervariable region (partial sequence) 7 , complete mtDNA control region 8,39 , mitogenome 8,31 , and whole-genome data 11,20,22 provided important insights in resolving the chicken phylogeny. In addition, recent genome-wide phylogenetic inferences provided a new perspective of wild species ancestry (i.e., G. g. spadiceus) of domestic chickens in southwestern China and Southeast Asia 11 . However, topological discrepancies have also been documented in genome-wide data, often explained by differences in data sources and taxon sampling 20,22,46 . The scope of the present study defines new evidence for modern chicken genetic information with increased data sources spanning Southeast Asia and Oceania. Furthermore, zooarchaeological DNA analysis can further clarify the evolutionary history of chickens in this region 12,47 .
Population genetic and phylogenetic analyses of more than 500 complete mtDNA control region sequences unveiled new perspectives on the population dynamics of SEA and Pacific chickens. Consistent with reports from various population genetic analyses, haplogroups A and B were widely distributed in East and Southeast Asia, while haplogroup E had the widest global distribution 7,8,19,31,40 . Haplogroup F was primarily represented in Myanmar chickens and shared this matriline with chicken populations in adjacent Yunnan Province, China 7,31,34 . Consistent with the phylogenetic analyses, the pairwise F ST value of Myanmar chickens was not genetically different from those of Yunnan chicken populations (Fig. 3a-b; Supplementary Table S6). This can be explained by the geographic proximity and the course of the Burma Road, which connects Myanmar and Yunnan Province 48 . Genetic differentiation of populations and PCoA analyses revealed genetic substructure between geographically isolated populations, i.e., between MSEA and ISEA chickens, South Asian and East Asian chickens, and South Asian and ISEA chickens ( Fig. 3a; Supplementary Table S6). Transregional population substructure was also observed within Southeast Asian chickens, reflecting deep phylogeographic diversification. Strong topological supports consistently define major haplogroup nomenclatures and provide evidence for the presence of a haplogroup D ancestral lineage (i.e., sub-haplogroup D2) from MSEA populations. A new matrilineage (i.e., sub-haplogroup V2) gave rise to the population of domestic chickens sampled in Cambodia, Laos, and Thailand, whereas their ancestral lineage (i.e., sub-haplogroup V1) was represented in Thai red junglefowl (i.e., G. g. gallus). The previously reconstructed mtDNA phylogenetic tree described by Huang et al. 31 assigned some of the previously identified haplogroup C samples to haplogroup V and linked them as a sister clade to the macrohaplogroup CD. However, because of the expanded sample distribution and increase in samples, we characterized haplogroup V as a sister group to haplogroup C only (Figs. 1, 2; Supplementary Fig. S2). This resulted in a clearer  www.nature.com/scientificreports/ reclassification of macrohaplogroup CDV (Fig. 1b). Interestingly, the ancestral matrilines classified under subhaplogroup D2 and haplogroup V were identified in sampling areas along the Lower Mekong subregion, for example, in Champasak and Bolikhamsai provinces in Laos, in Kampong Cham, Mondulkiri, Stung Treng, and   Table S1). The favorable climatic conditions and vegetation in this area are suitable for the red junglefowl (and their earlier descendants) to diversify and expand their distribution within their native range 6,10,49 . Collias & Saichuae 50 observed that the bird is drawn to primitive agriculture and converted primary forest. The bird also thrives and populates in the bamboo forest with lower elevations as well as near water holes or streams. In addition, migratory junglefowls have been sighted in the areas closer to the Mekong River, apparently attempting to cross it: "In crossing, the birds fly up as high as they can go, and then attempt to glide across… This movement does not seem to be caused by lack of food as the birds are extraordinary plump and in good condition. It is not easy to understand why it is taking place, as conditions on both sides of the Mekong seem the same" 51 . As one of the most geologically dynamic regions in the world, the Indo-Burma Biodiversity Hotspot has the highest recorded bird species (> 1,200) in the entire Asia-Pacific region 14,17 . The favorable seasonal weather patterns (i.e., dry northeast monsoon) and vegetation in much of the south, central, and west of the hotspot 17 make it a suitable habitat for chicken dispersal 10,12,18 . Meanwhile, the presence of the sub-haplogroup D1b classified for the Philippine and Pacific chicken populations is well documented in the present study and strongly supported by bootstrap and posterior probabilities. This matriline represents genetic uniformity and shows no significant signals of population structure despite geographic isolation between the Philippines and the Pacific region 36,37 . This recently expanded lineage is unique to this region, suggesting a human-mediated scenario of its phylogeography. This may be due to the dispersal of Austronesian speakers to the Philippines (ca. 4000 cal. BP) and continued movement eastward to the Melanesian islands (ca. 3300-3150 cal. BP) and as far as Remote Oceania [52][53][54] . This translocation route has been reliably defined by the recovered ancient DNA from Polynesian chickens, which identified the Philippines as a homeland for the diversity of Pacific chickens 36,55 . Similarly, the phylogeographic dispersal of the sub-haplogroup D1b, which first diversified in the Philippine archipelago, likely corresponds to the initial introduction pattern of their ancestral matriline (i.e., sub-haplogroup D2) from MSEA. This translocation pattern may have been influenced by the numerous waves of human migration to the Philippines brought by the Negritos across the continental landmass of Sundaland [56][57][58] . The introduction of the Manobo and Sama ancestry into the southern Philippines and Palawan cannot be ruled out, as they showed high genetic relatedness to MSEA-affiliated populations 58 . The timing of migration of people of Manobo ancestry (> 12,000 years ago) and people of Sama ancestry (~ 8000-12,000 years ago) is the closest possible translocation scenario 58 , which is consistent with archaeological evidence suggesting that the domestication of chickens in Southeast Asia occurred long before 8000 BP 6 . However, there are few reports of chicken remains in Southeast Asia 24 , and prehistoric exploitation has yet to be discovered 25 . Therefore, zooarchaeological and paleoclimatic studies are essential to identify their exact geographic center of origin reliably. On the contrary, we cannot assume a unidirectional north-to-south translocation of chickens from Taiwan because Taiwan's indigenous chickens (e.g., Ju-Chi) and gamecock (Hua-Tung) share genetic similarities with East Asian chicken haplotypes and populations introduced from Japan and the Indian subcontinent 59 .
Our coalescent-based Bayesian demographic analyses detected earlier effective population size expansion in MSEA chickens, while island populations showed more recent demographic growth signatures. Although our BSP results consider relevant sampling schemes with high sample sizes per demes, we still carefully acknowledge the potential impact of population structure on demographic estimates 60 . The timing of the demographic expansion of MSEA chickens observed in this study can be explained by the cultural importance of stock-raising in the archaeological sites of Non Nok Tha and Ban Chiang in Thailand around ca. 4000-3000 BP 49 . Bones of animals (e.g., pig, cattle, dog, deer, and chicken) and clay animal figurines were excavated in the human burial sites, suggesting that animals were part of the ritual practices during prehistoric inhumation 49 . It was well documented that agriculture and animal-raising were among the subsistence activities of domestic communities during prehistoric settlements in the broad valleys of the Lower Mekong 49 . In addition, ancient DNA of Thai chickens that were recovered in Ban Non Wat dated back to around 2500 BP, also supported the demographic expansion of MSEA chickens 24 . Recent morphological bone identification further documented the existence of chicken remains from other known archaeological sites in Thailand as early as 4000 BP 25 . On the other hand, the demographic expansion pattern of the island chicken population seems to suggest the timeline of Austronesian settlement in the region 61,62 .
The time tree phylogeny in the coalescent framework allowed us to estimate nodal ages of haplogroups relevant to this study. We combined primary calibration (i.e., fossil) from ancient Pacific chickens 36 and secondary calibration from previous estimations by Lawal et al. 20 . The latter calibration can provide close derived estimates from true time depending on the type of primary calibrations used 63,64 . Modelling the minimum-maximum constraints allows proximate measurement of uncertainties for estimated times and includes true time boundaries in the derived time CI range [63][64][65][66] . These calibrations and our robust phylogenetic trees allowed us to estimate the divergence of major haplogroups and the coalescence ages of some lineage-specific matrilines that shaped the population demographics of Southeast Asian chickens. For example, the coalescence time estimate for the node of macrohaplogroup CDV is estimated around 6.67 kya (95% HPD: 4235-7996 years). A similar age estimate was reported by Huang et al. 31 under a relaxed molecular clock model using the same molecular rate. The evidence of earlier coalescence age of haplogroup D ancestral matriline (i.e., sub-haplogroup D2) from MSEA exemplified dispersal patterns to the ISEA, and thereafter island clade diversified as a distinct group, a phylogeographic scenario that was also documented in other avian taxa 15,67,68 . Earlier paleoenvironment and biogeographic evidence 69 and more recent evidence on stable carbon isotope records from bat guano sequences 70  www.nature.com/scientificreports/ coalescence age (3.9 kya; 95% HPD: 2125-5880 years) than the previous estimate of this reclassified haplogroup 31 . Caution is warranted for this interpretation because the coalescence age estimate of gene copies in ancestral populations is not equivalent to a population split 72,73 , nor does it represent the actual onset of domestication.
In conclusion, this study provides a comprehensive insight into the genetic diversity and unique population dynamics of Southeast Asian chickens. High-resolution matrilineal phylogeny sheds new light on the evolutionary history of globally acknowledged haplogroups of SEA and Pacific chickens. It provides evidence of a new divergent matrilineage (i.e., haplogroup V) distributed across its native range in the Lower Mekong subregion. The phylogeographic and time tree phylogeny suggests human-mediated translocation of the haplogroup D ancestral matriline (i.e., haplogroup D2) from MSEA, which later diversified, forming a divergent sub-haplogroup D1b distinct to the island populations (i.e., Philippine-Pacific subclade). Future integrated genome-wide and environmental adaptation studies are required to unravel new elements of genomic evolution of SEA chickens for sustainable genetic improvement for climate resilience, effective management strategies, and future conservation endeavors.  Supplementary Table S1. Sampling was carried out from the unrelated individuals (e.g., sampling from a different known family and different sites within every province) to avoid lineage contamination during the later analyses. Genomic DNA was extracted from stored whole blood samples using the phenol-chloroform method 74 .

Materials and methods
The final dataset was complemented with previously published sequences of Philippine chickens (n = 129) 37 and directly submitted sequences of Indonesian (n = 10) and Pacific chickens (n = 11) retrieved from GenBank (Supplementary Table S1). Sequence alignment. Three hundred sixty-nine complete mtDNA control region sequences generated in this study were initially edited using GeneStudio Pro tool (GeneStudio, Inc., http:// genes tudio. com/). Ambiguous sites were trimmed, and cleaned sequences were aligned in MEGAX 77 with ClustalW 78 . Aligned nucleotide sequences were viewed using BioEdit 7.2.5 software 79 . All newly generated sequences were deposited in the GenBank database with accession numbers OM240181-OM240549 (Supplementary Table S1).

Genetic diversity and phylogenetic inference. Intrapopulation level and intraclade genetic diversity
indices such as the number of haplotypes (Ht), haplotype diversity (Hd), and nucleotide diversity (π) were estimated using the DnaSP v6.0 software 80 .
Phylogenetic analyses were inferred using two different model-based approaches: maximum-likelihood (ML) and Bayesian inference (BI). Maximum-likelihood analysis was performed in IQ-TREE 81 with the bestfit substitution model, TIM2 + F + I + G4, based on the Bayesian Information Criterion (BIC) determined by Modelfinder 82 . Statistical node support was calculated using ultrafast bootstrap support 83 and SH-aLRT 84  Population genetic structure and demographic inference. The population pairwise net genetic distance based on population pairwise F ST (significant values were accepted at p < 0.05) was estimated using Arlequin v3.5.2.2 software (with 10,000 permutations) 89 . Population pairwise F ST values were plotted into the principal coordinate analysis (PCoA) using GenAlEx v6.503 90 to visualize the pattern of genetic relationships between geographical populations. Estimation of the genetic structures was calculated by the analysis of molecular variance (AMOVA) as implemented by Arlequin v3.5.2.2 software. The level of significance was evaluated based on 1,000 random permutations.
Inference for the population growth model was initially estimated by statistical neutrality tests, such as Tajima's D 91 and Fu's Fs statistics 92 . These population expansion tests measure haplotype frequencies under neutrality. Statistical tests and confidence intervals were based on a coalescent simulation algorithm under a neutral infinite-site model. To further support the inference for the population expansion signal, a coalescentbased Bayesian Skyline Plot (BSP) was cautiously used to quantify the relationship between genealogies and the demographic history of the population 93 . BSP was simulated to infer deeper insights into the demographic history of Southeast Asian and Pacific chicken populations as implemented in BEAST v2.6.3 85 . BSP was generated with a relaxed molecular clock model and setting with 3.13 ×10 −7 mutations/site/year rate 86 . The piecewise constant function and HKY + G4 nucleotide substitution model as determined by BIC in jModelTest v2.1.1 94 was used for the analysis. The MCMC chain was run for 5 ×10 7 generations, with a sampling of parameters every 5000 steps and 5 ×10 6 generations served as burn-in. Convergence of the posterior estimates of the effective population size (N e ) to the likelihood stationary distribution was evaluated using Tracer v1.7.1 software 95 .
Divergence time estimate. Bayesian analyses were performed to estimate divergence times using the program BEAST2 v2.6.6. We employed a relaxed molecular clock model, which allows substitution rates to vary across branches setting with 3.13 × 10 −7 mutations/site/year rate 86 under uncorrelated lognormal distribution and GTR + G4 substitution model as determined by BIC in jModelTest v2.1.1. We set a coalescent-based constant population to model the tree prior. The ancient DNA records of Polynesian chickens were used to calibrate the crown node of sub-haplogroup D1b (Philippine-Pacific sub-clade) (Supplementary Table S1). For this calibration point, we used a lognormal prior (mean: 2.5, SD: 0.20, offset: 0) with the maximum age of the archaeological record set as the minimum bound for the crown calibration 63,65,66,96 . For the calibration of the root node of the tree, we used the established divergence time between red junglefowl and domestic chickens (8093 years CI: 7014-8768) 20 as a secondary calibration. We used a lognormal prior (mean: 8.09, SD: 0.05, offset: 0) covering the confidence interval range of the divergence time estimate 63,66 . Time tree analysis was run for 50 million generations, sampling every 5000 generations, and the initial 10% trees of each MCMC run were discarded as burn-in. The resulting log files were examined in Tracer v1.7.1 software 95 to confirm acceptable mixing and convergence of all parameters in the independent runs and adequate effective sample sizes (ESS > 200). The MCC tree was created from the tree file using TreeAnnotator v2.6.3 85 with the posterior probability set to 0.5 and common ancestor node heights summarized. These results were visualized as a single tree in FigTree v1.4.4. (http:// tree. bio. ed. ac. uk/ softw are/ figtr ee/).

Data availability
The complete mtDNA D-loop sequences are deposited and available in GenBank database (Accession Numbers: OM240181-OM240549). www.nature.com/scientificreports/